# -*- coding: utf-8 -*-
"""
Created on Mon Jan 18 23:27:02 2021

@author: Lenovo
"""

import numpy as np
from numpy import random
import matplotlib.pyplot as plt
global Gm,F0,Np,CR,D   #global是全局变量的标识符
import math

# 加载数据集
dataMat=np.loadtxt(open('D:\\冠心病稳定期VD.csv'),delimiter=',',skiprows = 1)
row=dataMat.shape[0] #数据集行数，数据条数
col=dataMat.shape[1] #数据集列数，症状个数
print(row,col)
# 根据选择生成相应约束
E1=0.5#主症的阈值
E2=0.1#次症的阈值
###################DE参数####################
Gm=200#进化代数
G=0#当前代数
D=2*col+1#参数维度
Np=8*D#种群规模
μ=Np
λ=3*μ
CR=0.2#交叉概率
F0=0.5#变异概率
ξ=200#自定义约束标准差值
tk=21#约束条件数目
pm = 0.05
mutation=np.zeros(D)#变异得到的向量 np.zeros方式生成的是二维形式
u=np.zeros(D)#交叉得到的向量注意第二种变异策略不需要交叉
value_fitness=np.zeros(μ+λ)#第二选择策略的最终适应值
value_fitness_index=0
bestX = np.zeros(D) # 种群中最好的个体
bestX_list=[]
#################函数定义与声明####################
def PX(v):#约束条件
    px=[]
    px.append(E1-((v[0]-v[1])/(1-v[1])))#5个主症，满足每个主症的公式大于阈值
    px.append(E1-((v[18]-v[19])/(1-v[19])))
    px.append(E1-((v[20]-v[21])/(1-v[21])))
    px.append(E1-((v[22]-v[23])/(1-v[23])))
    px.append(E1-((v[24]-v[25])/(1-v[25])))
    px.append(((v[2]-v[3])/(1-v[3]))-E2)#8个次症，满足每个主症的公式小于阈值
    px.append(((v[4]-v[5])/(1-v[5]))-E2)
    px.append(((v[6]-v[7])/(1-v[7]))-E2)
    px.append(((v[8]-v[9])/(1-v[9]))-E2)
    px.append(((v[10]-v[11])/(1-v[11]))-E2)
    px.append(((v[12]-v[13])/(1-v[13]))-E2)
    px.append(((v[14]-v[15])/(1-v[15]))-E2)
    px.append(((v[16]-v[17])/(1-v[17]))-E2)
    px.append(-((v[2]-v[3])/(1-v[3])))
    px.append(-((v[4]-v[5])/(1-v[5])))
    px.append(-((v[6]-v[7])/(1-v[7])))
    px.append(-((v[8]-v[9])/(1-v[9])))
    px.append(-((v[10]-v[11])/(1-v[11])))
    px.append(-((v[12]-v[13])/(1-v[13])))
    px.append(-((v[14]-v[15])/(1-v[15])))
    px.append(-((v[16]-v[17])/(1-v[17])))
    return px

def F1(v):#目标函数
    allsumlikelood=0#个体的BIC分,等于所有条数据的似然度之和再减去罚分
    sumlikelood=[]#个体在每条数据上的似然度
    #penalize是计算BIC分需要的后半部分,参数个数D除以2然后乘以log(row)
    penalize=(D/2)*(math.log(row))#对模型复杂程度的罚分
    for i in range(row):#对数据集中的每条数据
        #先算p(z=1)情况下概率
        likeloodforZis1=v[D-1]#p(z=1)
        likeloodforZis0=(1-v[D-1])#p(z=0)
        for j in range(col):
            if dataMat[i][j]==1:#p(x=1|z=1)时
                likeloodforZis1=likeloodforZis1*v[2*j]
                likeloodforZis0=likeloodforZis0*v[2*j+1]
            else:#p(x=0|z=1)时
                likeloodforZis1=likeloodforZis1*(1-v[2*j])#p(x=0|z=1)
                likeloodforZis0=likeloodforZis0*(1-v[2*j+1])#p(x=0|z=0)
        sumlikelood.append(likeloodforZis0+likeloodforZis1)
    allsumlikelood=sum(np.log(sumlikelood))
    return (allsumlikelood-penalize)

def G1(v):#约束函数1
    px=PX(v)
    g=[max(0,con) for con in px]
    gm=[pow((1.0*g[i]/p_dim[i]),2) for i in range(len(g))]
    gf=sum(gm)*1.0/tk
    return gf

def G2(v):#约束函数2
    px=PX(v)
    g=sum([pow(max(0,con),2) for con in px])
    return g

def getbest():#得到种群中最优的个体
    if(ψ==0):#全为不可行解
        return XG[np.argmin(Vcon)]
    elif(ψ==1):#全为可行解
        return XG[np.argmax(Vfun)]
    else:#可行解与不可行解皆有
        maxp=Vfun[np.argmin(Vcon)]#计算最优函数值
        return XG[Vfun.index(maxp)]#返回最优位置

def randn(i,number):#生成随机数
    randset=set()#创建集合，无序不重复
    if(number==3):
        while len(randset)!=4:
            randset.clear()  # 清空set
            randset.add(i)
            # 随机数范围为[0,size-1]的整数，三次
            r1=random.randint(0,Np-1)
            r2=random.randint(0,Np-1)
            r3=random.randint(0,Np-1)
            randset.add(r1)  # 添加进set，保证每个数都不相同
            randset.add(r2)
            randset.add(r3)
        return r1,r2,r3
    elif(number==5):
        while len(randset)!=6:
            randset.clear()  # 清空set
            randset.add(i)
            r1=random.randint(0,Np-1)
            r2=random.randint(0,Np-1)
            r3=random.randint(0,Np-1)
            r4=random.randint(0,Np-1)
            r5=random.randint(0,Np-1)
            randset.add(r1)
            randset.add(r2)
            randset.add(r3)
            randset.add(r4)
            randset.add(r5)
        return r1,r2,r3,r4,r5

def BGA(a,b,j): # 计算BGA变异所需变量
    rang=(a-b)*((1-G/Gm)**6)#a为上界，b为下界
    sum=0
    for k in range(D):
        rand=random.random()
        α=1 if rand<(1/D) else 0
        sum+=α*(2**(-k))
    return rang*sum

def Fitness(x):
    global value_fitness_index
    value=value_fitness[value_fitness_index]
    value_fitness_index+=1
    return value

###############种群初始化################
#---------佳点集方法---------
max0=1
min0=0
X0 =random.random(size=(Np,D))
for i in range(D):
    X0[:,i]=X0[:,i]*(max0-min0)+min0
XG=X0
#----------约束程度标准化-----------
P0=[]
for v in XG:
    P0.append(PX(v)) #将一维的px加入到二维的P0
P0=np.array(P0)#约束惩罚度总表
p_dim=np.max(P0,axis=0)
δ=np.max(p_dim)-np.min(p_dim)#标准化之后的差值
Vfun=[F1(v) for v in XG]#目标函数值Fx
if δ>ξ:#选择约束函数计算方式Gx
    Vcon=[G1(v) for v in XG]
else:
    Vcon=[G2(v) for v in XG]
#-------------准备工作--------------
num_1=μ-np.count_nonzero(Vcon)
ψ=num_1/μ#种群中可行解占的比例,代代更新
bestX=getbest()
#####################开始进化#################
while G<Gm:
    print(G,'代')
    G=G+1##最后移到后面    
    suanzi=math.exp(1-float(Gm)/(Gm+1-G))
    F=F0*(2**suanzi)
    Qt=[]
    #下面执行变异和交叉
    for i in range(Np):
    #--------------第一种变异策略 rand/1 ----------------------
        r1,r2,r3=randn(i,3)
        mutation=XG[r1]+F*(XG[r2]-XG[r3])
        for j in range(D):  # 判断变异后的值是否满足边界条件，不满足需要重新界定
            if mutation[j]>max0:
                mutation[j]=2*max0-mutation[j]
            elif mutation[j]<min0:
                mutation[j]=2*min0-mutation[j]
            else:
                mutation[j]=mutation[j]
        # 交叉
        randx=list(range(0,D-1,2))
        random.shuffle(randx)
        for j in range(0,D-1,2):
            if random.rand()>CR and randx[0]!=j:
                u[j]=XG[i,j]
                u[j+1]=XG[i,j+1]
            else:
                u[j]=mutation[j]
                u[j+1]=mutation[j+1]
        if np.random.rand()>CR:
            u[D-1]=XG[i,D-1]
        else:
            u[D-1]=mutation[D-1]
        Qt.append(u)
    #--------------第二种变异策 改进的current-to-best/ ---------------
        u = np.zeros(D)  # 对变量u重新赋值初始化
        rand=random.rand()
        r1,r2,r3=randn(i,3) # 此处重新生成r1,r2,r3
        if(rand<ψ):#rand()随机生成一个处于0~1的数
            mutation=XG[i]+rand*(XG[r1]-XG[i])+F*(XG[r2]-XG[r3])
        else:
            mutation=XG[i]+rand*(bestX-XG[i])+F*(XG[r1]-XG[r2])
        if(rand<pm): # 下面执行改进的BGA变异
            rand1=random.rand() # 重新生成rand随机数
            if(rand1<(1/D)):
                for j in range(D):
                    rand2=random.rand()
                    bga=BGA(max0,min0,j)
                    mutation[j]=mutation[j]+bga if rand2>0.5 else mutation[j]-bga
        for j in range(D):  # 判断变异后的值是否满足边界条件，不满足需要重新界定
            if mutation[j]>max0:
                mutation[j]=2*max0-mutation[j]
            elif mutation[j]<min0:
                mutation[j]=2*min0-mutation[j]
            else:
                mutation[j]=mutation[j]
        u=mutation
        Qt.append(u)
    #--------------  第三种变异策略 rand/2变异--------------
        u = np.zeros(D)
        r1,r2,r3,r4,r5=randn(i,5)
        mutation=XG[r1]+F*(XG[r2]-XG[r3])+F*(XG[r4]-XG[r5])
        for j in range(D):  # 判断变异后的值是否满足边界条件，不满足需要重新界定
            if mutation[j]>max0:
                mutation[j]=2*max0-mutation[j]
            elif mutation[j]<min0:
                mutation[j]=2*min0-mutation[j]
            else:
                mutation[j]=mutation[j]
        # 交叉
        randx=list(range(0,D-1,2))
        random.shuffle(randx)
        for j in range(0,D-1,2):
            if random.rand()>CR and randx[0]!=j:
                u[j]=XG[i,j]
                u[j+1]=XG[i,j+1]
            else:
                u[j]=mutation[j]
                u[j+1]=mutation[j+1]
        if np.random.rand()>CR:
            u[D-1]=XG[i,D-1]
        else:
            u[D-1]=mutation[D-1]
        Qt.append(u)
    #------------下面进行选择操作---------------
    QCV=[G2(v) for v in Qt] # 变异种群Qt的违约值
    num_2=num_1+(λ-np.count_nonzero(QCV)) # 混合种群中的可行解数量
    fp=num_2/(μ+λ)
    #print('num_1=',num_1,'ψ=',ψ)
    print('num_2=',num_2,'fp=',fp)
    hy=np.vstack((XG,np.array(Qt))) # 形成混合群体
    # 以下为不可行种群选择策略-----------
    sorted_hy=np.array(sorted(hy,key=G2)) # 根据违约值排过序的混合种群
    if(fp==0):
        P_temp=sorted_hy[:int((μ + λ)/2)]  # 对混合群体进行排序 择优录取原先的一半(μ+λ)/2
        half_P_temp=P_temp[:int(μ/2)]  # 从P_temp中取得的μ/2个优秀个体
        remain_P_temp=P_temp[int(μ/2):]  # P_temp剩下的个体 需要随机从中取得μ/2
        rand_li=list(range(int(λ/2)))
        random.shuffle(rand_li)  # 构造随机数列
        XG=np.vstack((half_P_temp,remain_P_temp[rand_li[:int(μ/2)]]))  # 两部分合并组成新种群
    # 以下为可行解不可行解混合种群选择策略
    elif(0<fp<1):
        Z1=sorted_hy[:num_2]  # 可行种群，记为Z1
        Z2=sorted_hy[num_2:]  # 不可行种群
        value_fun_1=[F1(v) for v in Z1]#f(Z1)
        X1_best_value=np.min(value_fun_1)#可行种群中最优解与最差解
        X1_worst_value=np.max(value_fun_1)  
        value_con_2=[G2(v) for v in Z2]  # 不可行种群的违约值g(Z2)
        min_value_con_2=np.min(value_con_2)  # 不可行种群违约值的最大与最小
        max_value_con_2=np.max(value_con_2)
        value_fun=np.zeros(μ+λ)#f(Z1+Z2)'
        fnor_value=np.zeros(μ+λ)#f nor改进的目标函数值
        gnor_value=np.zeros(μ+λ)#g nor改进的违约值
        for i in range(sorted_hy.shape[0]):
            if i < num_2:
                value_fun[i]=value_fun_1[i]
            else:
                value_fun[i]=max(ψ*X1_best_value+(1-ψ)*X1_worst_value,F1(sorted_hy[i]))
                gnor_value[i]=(value_con_2[i - num_2]-min_value_con_2)/(max_value_con_2-min_value_con_2)  # Z2时的计算，Z1为0，可省略。
        fnor_value=(value_fun-np.min(value_fun))/(np.max(value_fun)-np.min(value_fun))# 改进目标函数值的标准化
        value_fitness=fnor_value+gnor_value# 最终的适应函数值
        sort_by_fit_hy=np.array(sorted(sorted_hy,key=Fitness))  #对混合种群根据适应值排序
        XG = sort_by_fit_hy[λ:]
    # 以下为可行解模式
    elif(fp==1):
        value_func=np.array(sorted(hy,key=F1))  #对混合种群根据目标函数值排序
        XG=value_func[λ:]
    # 本次进化结束，重新计算相关参数
    Vfun=[F1(v) for v in XG] #目标函数值Fx
    Vcon=[G2(v) for v in XG]  #种群的违约值
    num_1=μ-np.count_nonzero(Vcon)
    ψ=num_1/μ
    bestX=getbest()
    print(bestX)
    print('解',F1(bestX))
    print('约束',G2(bestX))
    bestX_list.append(F1(bestX))
    value_fitness_index=0
    #进化结束
with open('(μ+λ)-DE.txt','w') as f:
    for i in bestX_list:
        print(i, file=f)
plt.plot(bestX_list)
plt.show()
print("最优解：%f" %F1(bestX))
print("最优位置：",bestX)
print("违约值：%f" %G2(bestX)) 

best_vector=bestX
thefirst=math.log((1-best_vector[D-1])/best_vector[D-1])#计算阈值前半部分
sumv=0#计算阈值公式中的后半部分
for i in range(0,D-1,2):
    temp=(1-best_vector[i])/(1-best_vector[i+1])
    sumv=sumv+math.log(temp)
thesorld=thefirst-sumv
scalerate=thesorld/10#缩放率等于原始阈值除以10，这样把阈值统一为10
###########下面计算每个症状的分值
name=['胸闷','肢体困重','体胖','口黏','面色浑浊','大便黏滞','脘腹痞满','嗜睡','纳差','舌胖边有齿痕','苔腻','苔白滑','脉濡或滑']
score = 0
for i in range(0, D - 1, 2):
    s=(best_vector[i]/(1-best_vector[i]))/(best_vector[i + 1]/(1-best_vector[i+1]))
    score=(math.log(s))/scalerate#原始分值除以缩放率，这样把阈值统一为10
    j=(int)(i/2)#症状有13个，变换为症状的下标
    print('%s:'%(name[j]),end='')
    print('%0.2f'%score)#小数点后保留2位
